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Abstract: On directed Barabasi- Albert and Small- World networks the Ising 
model with spin 5 = 1,3/2 and 2 is now studied through Monte Carlo simulations. 
In this model, the order-disorder phase transition of the order parameter is well 
defined on Small- World networks for Ising model with spin S = 1. We calculate the 
value of the critical temperature Tc for several values of rewiring probability p of 
the directed Small- World network. This model on directed Small- World networks 
we obtained a second-order phase transition for p = 0.2 and first-order phase 
transition for p = 0.8. The critical exponentes P/i', ^jv and Xjv were calculated 
for p = 0.2. On directed Barabasi- Albert we show that no there is phase transition 
for Ising model with spin S = 1, 3/2 and 2. 
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Introduction 

This paper deals with Ising spin on directed Barabasi- Albert (BA) and 
Small-World(SW) networks. Sumour and Shabat [U [2] investigated Ising 
models with spin S = 1/2 on directed B A networks [3] with the usual Glauber 
dynamics. No spontaneous magnetisation was found, in contrast to the case 
of undirected BA networks [H [5], [6] where a spontaneous magnetisation was 
found below a critical temperature which increases logarithmically with sys- 
tem size. In S=l/2 systems on undirected, scale-free hierarchical-lattice SW 
networks [7] , conventional and algebraic (Berezinskii-Kosterlitz-Thouless) or- 
dering, with finite transition temperatures, have been found. Lima and Stauf- 
fer [Sj simulated directed square, cubic and hypercubic lattices in two to five 
dimensions with heat bath dynamics in order to separate the network effects 
form the effects of directedness. They also compared different spin flip al- 
gorithms, including cluster flips [H], for Ising-BA networks. They found a 
freezing-in of the magnetisation similar to [H |2] , following an Arrhenius law 
at least in low dimensions. This lack of a spontaneous magnetisation (in the 
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usual sense) is consistent with the fact that if on a directed lattice a spin Sj 
influences spin Si, then spin Si in turn does not influence Sj, and there may 
be no well-defined total energy. Thus, they show that for the same scale-free 
networks, different algorithms give different results. The g-state Potts model 
has been studied in scale-free networks by Igloi and Turban [10] and depend- 
ing on the value of q and the degree-exponent 7 first- and second-order phase 
transitions are found, and also by Lima [TT] on directed BA network, where 
only first-order phase transitions have being obtained independent of values 
of q for values of connectivity z = 2 and z = 7 of the directed BA network. 
More recently, Lima [12] simulated the the Ising model for spin S = 1 on 
directed BA network and different from the Ising model for spin S = 1/2, 
the unsual order-disorder phase transition of order parameter was seen; this 
effect is re-evaluated in the light of the time depemdence presented below. 
We study the Ising model for spin 1, 3/2 and 2 on directed BA network. The 
Ising model with spin 1, 3/2 and 2 was seen not to show a usua/ spontaneous 
magnetisation and this decay time for flipping of the magnetisation followed 
an Arrhenius law for HeatBath algorithms that agree with the results of the 
Ising model for spin S' = l/2[Tl[2]on directed BA network. Edina and Lima 
[T3j calculate the exponentes P/p, ^jv, and 1/z/ exponents for majority- vote 
model on directed SW networks of Sanchez et al. [H] , and on these networks 
the exponents are different from the Ising model a two-dimensional and in- 
depend on the values of rewiring probability p of the directed SW networks. 
Here we also study the Ising model for spin = 1 on directed SW, we ob- 
tained a second-order phase transition for p = 0.2 and a first-order phase 
transition for p = 0.8. We also calculate the critical exponents iijv, 7/1/ and 
1/z/ for p = 0.2. 



Model and Simulator! 

Ising model on directed BA Networks: 

We consider the spin 1,3/2 and 2 Ising models on directed BA Networks, 
defined by a set of spins variables 5* taking the values ±1 and for 5=1, 
±3/2 and 1/2 for 5* = 3/2, and ±1, ±2 and for 2, respectivaly, located on 
every site of a directed BA Networks with sites. 
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Figure 1: Reciprocal logarithm of the relaxation times on directed BA net- 
works for 5* = 1 to 5" = 2. 

The probability for spin Si change its state in this directed networks is 
Pi = 1/1 + eM'^Ei/kBT), Ei = -jY^SiSk (1) 

k 

where A;- sum runs over all nearest neighbors of S;,. In this network, each new 
site added to the network selects with connectivity z = 2 already existing 
sites as neighbours influencing it; the newly added spin does not influence 
these neighbours. 

To study the spin 1,3/2 and 2 Ising models we start with all spins up, a 
number of spins equal io N — 500000, and Monte Carlo step (MCS) time up 
2, 000, 000, in our simulations, one MCS is accomplished after all the spins 
are updated, here, with HeatBath Monte Carlo algorithm. Then we vary the 
temperature and study nine samples. The temperature is measured in units 
of critical temperature of the square-lattice Ising model. We determine the 
time r after which the magnetisation has flipped its sign for the first time, 
and then take the median value of our nine samples. So we get different 
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Decay time tau for 7 million with m=2 (left +) and 2 million witfi m=7 (rigtit +); and smaller systems 
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Figure 2: Characteristic time for M(r) = 3/4 using 7 million spins for m = 2 
neighbours and 2 million spins for m = 7 neighbours (+). Ten, hundred, 
and thousand times smaller systems are denoted by x, star, and squares. We 
plot the median over nine samples in this log- log plot. The two straight lines 
have negative slopes 8 (left) and 12 (right). From 

values Ti for different temperatures. 

In the study the critical behavior this Ising model (with spins 1, 3/2 and 
2) we define the variable m = Y^iLi Si/N as normalized magnetisation. 

Our BA simulations, using the HeatBath algorithm, indicate that the 
spins = 1, 3/2 and 2 Ising model do not display a phase transition and the 
plot of the time 1/ In r versus temperature in Fig. 1 shows that our BA results 
for all spins agree with the modified Arrhenius law for relaxation time, defined 
as the first time when the sign of the magnetisation flips:!/ ln(r) oc T + .... 
This result agrees with the results of Sumour et al. for spin S = 1/2 [2J with 
7 million of sites, see Fig. 2 (figure retired of reference |2j|) and our results 
are more reliable than both Lima |12j and Sumour et al. because that 
this are longer simulation times. 
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Magnetisation versus Temperature, L=128 
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Figure 3: Magnetisation as a function of the temperature T, for = 16384 
sites. A second-order phase transition for values of p = 0.2 and a first-order 
phase transition for p = 0.8 . 

Ising model on directed SW Networks: 

We consider the Ising model with spin S* = 1, on directed SW Networks, 
defined by a set of spins variables S taking the values —1,0 and +1, situated 
on every site of a directed SW Networks with sites. 

The probability for spin Si to change its state in these directed network 
is given eq. (1) where the sum runs over all nearest neighbors of S^. In this 
network, created for Sanchez et al. [H], we start from a two-dimensional 
square lattice consisting of sites linked to their four nearest neighbors by 
both outgoing and incoming links. Then, with probability p, we reconnect 
nearest-neighbor outgoing links to a different site chosen at random. After 
repeating this process for every link, we are left with a network with a density 
p of SW directed links. Therefore, with this procedure every site will have 
exactly four outgoing links and varying (random) number of incoming links. 
To study the critical behavior of the model we use the HeatBath algorithm 
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Figure 4: Binder's fourth-order cumulant as a function of K. We have 
p = 0.2 (part a) and p = 0.8 (part b) for L = 8, 16, 32, 64 and 128. 



and define the variable e = E/N, where E = 2J2iEi and m = J2i=i Si/N. 
From the fluctuations of the e measurements we can compute: the average 
e, the specific heat C and the fourth-order cumulant of e, 

u{K) = [<E >U/N, (2) 

C{K) = K^N[< > - < e >\,, (3) 

= [1 - ^^^]»- (4) 

where ii' = J/ksT, with J = 1, and ks is the Boltzmann constant. Similarly, 
we can derive from the magnetization measurements the average magnetiza- 
tion, the susceptibility, and the magnetic cumulants, 

m{K) = [< \m\ >]av, (5) 

xiK) = KN[< m^>-< \m\ >X,, (6) 

U,iK) = [1 - (7) 
3 < \m\ >^ 

where < ... > stands for a thermodynamic average and [...]av square brackets 
for a averages over the 20 realizations. 

In order to verify the order of the transition this model, we apply finite- 
size scaling (FSS) [15]. Intially we search for the minima of e fourth-order 
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Figure 5: 2/3 — -Bj^mm for p = 0.2 (circles) and 0.8 (squares). 



parameter of eq. (4). This quantity gives a qualitative as well as a quanti- 
tative description of the order of the transition [16]. It is known [TTj that 
this parameter takes a minimum value -Bj,m,m at effective transition temper- 
ature Tc{N). One can show [H] that for a second-order transition limAr^oo 
(2/3 — Bi min) = 0, even at T^, while at a first-order transition the same limit 
measurng the same quantity is small and (2/3 — Bi min) 7^ 0. 

A more quantitative analysis can be carried out through the FSS of the 
C fluctuation Cmax, the susceptibility maxima Xmax and the minima of the 
Binder parameter Bi^min- If the hypothesis of a first-order phase transition is 
correct, we should then expect, for large systems sizes, an asymptotics FSS 
behavior of the form [T9l 



= ac + hcN + ... (8) 
= ^x + KN + - (9) 

.^,m^n^aB,+hB,N + ... (10) 

Therefore, if the hypothesis of a second-order phase transition is correct, 
we should then expect, for large systems sizes, an asymptotics FSS behavior 
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Figure 6: Binder's fourth-order cumulant of the magnetisation versus tem- 
perature for p — 0.2. 

of the form 

C = ae, + L"/Vc(x)[l + ...], (11) 

[<|m|>U = L-^/Vm(^)[l + ...], (12) 

X = L^/Vx(^)[l + -], (13) 

^ = lVv^(x)[1 + ...], (14) 

where Creg is a regular background term, u, a, /3, and 7, are the usual critical 
exponents, and fi{x) are FSS functions with 

x^{K- K^)L^I'' (15) 

being the scaling variable, and the brackets [1 + ...] indicate corretions-to- 
scaling terms. Therefore, from the size dependence of M and x we obtained 
the exponents [3/u and 7///, respectively. The maximum value of susceptibil- 
ity also scales as L^l^ . Moreover, the value of T for which x has a maximum. 
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Figure 7: In M{Tc) versus In L for p = 0.2. 

j'Xmax — Tc{L), is expected to scale with the system size as 

T,{L)^n + bL-^/\ (16) 

were the constant b is close to unity. Therefore, the relations (14) and (16) 
are used to determine the exponent l/u. 

We have performed Monte Carlo simulation on directed SW network with 
various values of probability P- For a given p, we used systems of size L = 8, 
16, 32, 64, and 128. We waited 50, 000 Monte Carlo steps (MCS) to make the 
system reach the steady state, and the time averages were estimated from 
the next 50, 000 MCS. In our simulations, one MCS is accomphshed after all 
the spins are updated. For all sets of parameters, we have generated 20 
distinct networks, and have simulated 20 independent runs for each distinct 
network. 

Results and Discussion 
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Figure 8: Plot of In ;^"*"^(L) (square) and Inx(Tc) (circle) versus In L for 
p = 0.2. 

In Fig. 3 we show the dependence of the magnetisation M on the temper- 
ature, obtained from simulations on directed SW network with L = 128 x 128 
sites and two values of probability p: p = 0.2 for second-order transition and 
p = 0.8 for first-order transition. In Fig. 4 we plot Binder's fourth-order 
cumulant of e for different values of L and two different values of p: p = 0.2 
in part (a) and p — 0.8 in part (b). In Fig. 5 the difference 2/3 — Bi min is 
shown as a function of parameter 1/A^ for probability p = 0.2 (circles) and 
= 0.8 (squares) obtained from the data of Fig. 4. In Fig. 6 we show Binder's 
fourth-order cumulant of the magnetisation versus temperature for p = 0.2. 
The temperature obtained is Tc — 1.890(6). Figs. 7 we plot the dependence 
of the magnetisation at T = Tc versus the system size L. The slopes of curves 
correspond to the exponent ratio (3/v oi according to Eq. (12). The results 
show that the exponent ratio j3/v at is 0.379(21). 

In Fig. 8 we display the scalings for susceptibility at T = Tc(L) , x(Tc(L)) 
(circles), and for its maximum amplitude, xIl""^ (squares), and the scalings 
for susceptibility at T = Tc obtained from Binder's cumulant, x(^c) versus L 
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Figure 9: Reciprocal logarithm of the relaxation times on directed SW 
networks for versus T for p=0.2(+) and 0.8(x). 



for probabihty p = 0.2. The exponents ratio '-f/u are obtained from the slopes 
of the straight lines. The exponents 7/1/ of the two estimates agree (within 
errors). The values obtained are ^/v — 1.252(13) (circles) and 1.264(34) 
(squares), respectively. Therefore we can use the Eq. (16) , for p = 0.2, 
obtain the critical exponent 1/z/, that is equal to 1.209(165) and 1.248(196) 
obtained of according to Eq. (14). To improve our results obtained above 
we start with all spins up, a number of spins equal to = 640000, and 
time up 2,000,000 (in units of Monte Carlo steps per spins). Then we vary 
temperature T and at each T study the time dependence for 9 samples. We 
determine the time r after which the magnetisation has flipped its sign for 
first time, and then take the median values of our nine samples. So we get 
different values Ti for different temperatures T. In Fig. 9 show that the 
decay time goes to infinity at some positive T value. This behavior ensures 
that there is a phase transition for Ising model spin 5" = 1 on directed SW 
network. 
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Conclusion 

In conclusion, we have presented the Ising model spin spins S = 1, 3/2 
and 2 on directed BA and directed SW network(S' = 1). The Ising model does 
not display a phase transition on directedBA for spins S = 1/2, 1, 3/2 and 2. 
In the directed SW network [13] this model presents a first- and second-order 
phase transition which occurs with probability p = 0.8 and 0.2, respectively. 
The exponents obtained for p = 0.2 are different from the exponents the 
Ising model on square lattice, that it suggests that these exponents belong 
to one another class of universality. 
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